Adsorption and Absorption Energies of Hydrogen with Palladium

Thermal recombinative desorption rates of HD on Pd(111) and Pd(332) are reported from transient kinetic experiments performed between 523 and 1023 K. A detailed kinetic model accurately describes the competition between recombination of surface-adsorbed hydrogen and deuterium atoms and their diffusion into the bulk. By fitting the model to observed rates, we derive the dissociative adsorption energies (E0, adsH2 = 0.98 eV; E0, adsD2 = 1.00 eV; E0, adsHD = 0.99 eV) as well as the classical dissociative binding energy ϵads = 1.02 ± 0.03 eV, which provides a benchmark for electronic structure theory. In a similar way, we obtain the classical energy required to move an H or D atom from the surface to the bulk (ϵsb = 0.46 ± 0.01 eV) and the isotope specific energies, E0, sbH = 0.41 eV and E0, sbD = 0.43 eV. Detailed insights into the process of transient bulk diffusion are obtained from kinetic Monte Carlo simulations.


Molecular Beam Flux Calibration
A procedure for the estimation of the number of molecules in a molecular beam pulse was developed using laser based multi-photon-ionization and spatial ion imaging. Comparable methods have been proposed previously. 1,2 The procedure aims to characterize the spatial molecule density produced by one pulse. The number of molecules in a volume element is given as ( , , ) d d d ( , , ) = ( , , )d = max d d d Eq. S1 In this notation is the molecule density (molecules/ m 3 ) at the peak of the molecular beam max density distribution. The shape of the distribution is characterized by the shape functions , and which are assumed independent of each other. The product of three Gaussians describes the beam S3 shape accurately. Each shape function is peak normalized, such that the density at the center equals . max In addition, an ionization efficiency function was defined in order to account for varying ( , , ) molecule detection probability along the laser propagation ( -axis). The laser ionization function returns the probability of ionization at a point in space where is the probability of ionization at the center of the laser focus. In the present experiment, max the laser ionization profile is narrow compared to the molecular beam pulse in and . Only the component, which is the propagation axis of the laser beam, is similar in width to the molecular beam.
The number of molecules , that are ionized by a laser pulse, is given by integration of Eq. S1 weighted by the ionization probability.
For example, if static gas is ionized by the laser instead of a molecular beam pulse, the shape functions , and are replaced by 1. The number of ionized molecules is then given by the laser ionization profile and the constant background pressure. However, if a molecular beam pulse is characterized, the number of ions in a volume element is proportional to both the laser ionization d d d probability and the molecule density in the volume element. We assume that the laser ionization probability is independent of molecular density.
The molecular density in the beam was determined through calibration of the ion detector with a static background gas of known density. Eq. S4 gives the ratio of ionized molecules in the molecular beam and in the calibration with the static gas

Eq. S4
The number of ionized molecules using static gas does not include , and as the density is BG constant in space. The absolute ionization probability of the laser cancels, as the same laser max intensity was used for molecular beam and static gas ionization. The laser focus is small compared to the molecular beam shape in directions orthogonal to the propagation axis ( and ), therefore ∫ with being the constant value which applies for the whole relevant integration range. Experimentally, the laser focus is placed in the center of the beam such that . For = = 1 this cannot be done as the laser focus and molecular beam are of comparable size along the -axis.
From Eq. S4 we obtain Eq. S5 To obtain the particle density of the molecular beam at its center , the density in the static gas MB max experiment needs to be known. It was obtained from the N 2 equivalent ion gauge reading ,

BG IG
which was converted to H 2 (D 2 ) pressure using the factor . 3 Hydrogen was H 2 (D 2 ) = 2.17 (2.86) supplied to the chamber by a leak valve. The established constant background pressure was used to determine from the ideal gas law where is Boltzmann's constant and is the absolute temperature in the chamber (296 K). The B ratio of produced ions when the laser hits the center of the molecular beam and when static gas is used is obtained as the experimental signal intensity ratio where is the MB / BG = MB / BG number of counts measured in arbitrary units.
The laser ionization function is obtained by comparison of the beam shape function and the ionized fraction of the molecular beam profile as = Eq. S7 The beam shape function is obtained from an image, where ion signal was integrated while the laser focus is translated along the -axis. This leads to a spatially homogeneous ionization efficiency, which reflects the molecular beam density shape on the ion image only. The quantity is the shape resulting from the laser focus positioned in the center of the molecular beam. We emphasize that only the shape of the quantities , and is of importance.
With all these quantities in hand, Eq. S5 can be evaluated to yield (molecules/m 3 ). In order to MB max obtain the number of molecules within one beam pulse Eq. S1 needs to be integrated over the entire space. Two more transformations are employed to evaluate the integral. i) The integration along (molecular beam propagation axis) is obtained by scanning the beam-laser delay time at fixed -S5 position of the laser focus. This leads to the temporal shape function . The mean velocity of the MB molecular beam is used to convert time to space. ii) The beam shape along cannot be measured MB in our detection geometry. However, radial symmetry can be assumed around the -axis, as the beam is produced from a round orifice and passes a round skimmer. Integration in the -plane is therefore performed after transformation into polar coordinates, taking as the radius. Finally, Eq. S8 results for the total number of molecules in a molecular beam pulse.
The number of molecules within a molecular beam pulse is distributed over the surface according to the molecular beam shape function . When the dimension of propagation or is integrated, the exposed amount of H 2 molecules is given as: where is the width of the molecular beam in or and is the density of ( = 1.53 × 10 15 The time-dependent adsorbing flux of H atoms resulting from the exposure of H 2 is then given by with as the dissociative sticking coefficient. This expression can be found in Eq. S1 of the 2 HH 0 HH 0 main text.

Coarse Grained Grid Diffusion
It is well known that hydrogen atoms easily penetrate into the bulk of Pd metal. Using the well-known concentrations decrease, such that it is no more necessary to simulate the concentration of every single layer. Instead, the concentration evolution is simulated in a grid with exponentially growing grid point separation with increasing depth in the crystal, see Figure S1. The correct description of the concentration exchange between a grid point and its neighbours is derived from the diffusion equation. In general, the 1D diffusion equation along the spatial axis is given as Eq. S11 where is the concentration and is the diffusion coefficient. Evaluation of the 2 nd derivative of the concentration at spatial point with concentration is approximated by finite differences using the left and right neighbouring points. 8 Eq. S12

S7
Further, the diffusion coefficient is replaced by its microscopic representation in one dimension: = .
Here, is the total jump frequency for a particle located in a discrete site. The site is connected Γ 2 /2 Γ to a left and a right neighbour by the jump length . is composed of two equivalent contributions Γ for hopping to the left and to the right . As describes the transition from one site to the next l r bulk in one direction, it follows: and . Finally, the position variable is expressed in Γ = 2 bulk = bulk 2 units of the jump length by substituting , where is the number of single layers needed to = × reach the depth at point . With , the diffusion equation is written as Eq. S13 In summary, Eq. S13 evaluates the rate of concentration change at grid point , where the spatial axis is given as depth below the surface in units of the single jump length. If grid points are separated by the (111) layer-width, the diffusion equation turns into the kinetic equation for explicit hopping in a uniform grid. Eq. S14 -Eq. S16 represent the full set of differential equations for H used to propagate the model in time. Analogously, equations for D are implemented. Typically, the grid consists of about 100 points (compare Figure S1).

Bulk Diffusion Constants for H and D
The bulk diffusion constant was obtained from an Arrhenius fit to experimental data from literature [5][6][7] (see Figure S2). The fit yields the diffusion barrier and the pre-factor .
Eq. S17 The microscopic hopping pre-factor is obtained from the macroscopic as Eq. S18 where is the distance between two (111) layers, which relates to the lattice constant of palladium pm 9 as . = 389 = 3 The resulting parameters are listed in Table S2:

Bulk Potential and Partition Function
The effective potential used in the desorption/diffusion model for an H atom absorbed in the bulk of Pd is shown in Figure S3 and is described by Eq. S19. bulk ( )  Eq. S19 Note that the potential function was obtained by fitting the absorption enthalpy data from previous work 10 together with the present kinetic traces as described in the main text. In the fitting procedure, we forced the harmonic frequency, calculated from the minimum energy position of the octahedral site (dashed line Figure S3), to the experimental value of 0.069 eV, that has been measured by neutron scattering. 11 This treatment guarantees an accurate ZPE and anchors the potential to measured quantities. The best fit curve obtained here is highly anharmonic and its properties are consistent with conclusions drawn from neutron scattering studies. 11 Further, we also incorporated a site that mimics the tetrahedral site. At high energy, the free particle limit is approached. This description features physical anchor points and is capable to effectively fit the experimental observations.
The partition function of is the most critical quantity to describe the temperature dependence bulk of the absorption enthalpy as well as being important to describe the surface to bulk transition rate.
It is obtained from the classical partition function corrected for low temperature quantum effects using the Pitzer-Gwinn correction. 12,13 The three-dimensional partition function is thus given as Eq. S20

Modelling the Absorption Enthalpy
The The absorption enthalpy is shown as a function of temperature in Figure 3 of the main paper.

Thermal Sticking Coefficient
Thermal sticking coefficients of H 2 and D 2 were obtained by averaging energy dependent sticking coefficients measured at the surface normal. 14,15 The applied procedure is described in Ref. 1 . Figure   S4 shows the result. The thermal sticking coefficient of HD was assumed to be the average of the values for H 2 and D 2 . We obtained a temperature independent value of . S HD 0 = 0.535 S12 Figure S4: Thermal sticking coefficient of H 2

H(D)-Recombination Rate Constant
The hydrogen atom recombination rate constant is given by Eq. S8 of the main text. The specific contributions to this formula are described in the following. For the hydrogen atom partition function (Eq. S23), separability of motion in the -plane from motion along the -coordinate is assumed. The in-plane partition function was obtained by solving the single particle Schödinger equation where is the ground-state energy and is chosen such that the partition function is converged for 0 a given temperature. The procedure is described in great detail in Ref. 1

. The in-plane H-Pd interaction
PES was obtained from DFT by calculating a representative -energy grid for each surface (see SI sec. 14). The -coordinate (out of plane) of the atom was held at the minimum energy at each point. This assumes that in-plane and out of plane movement are decoupled. Figure S5   . Therefore, we need to account for the coupling of z-stretch frequency to the in-plane coordinates. 16 The -dependent frequencies are used to obtain the mean partition PdA ( , ) function . This quantity is the spatial average of the vibrational partition function, weighted The partition function of the gas phase product molecule is given as In Eq. S26 denotes the molecular mass, is the symmetry number, is the rotational constant in the vibrational ground-state and is the vibrational frequency of the molecule.
Required molecular constants were taken from experiments and are reported in Table S4.
The partition functions were defined to start from the lowest quantum state and not from the classical minimum of the PES. Therefore, the binding energy includes the zero-point energies. It

Eq. S27
In Eq. S27 denotes the classical dissociative binding energy of hydrogen at the surface. Table S3 rec summarizes the zero-point energies of the adsorbed H-atom.  17 . These are used for the uncertainty evaluation.

Uncertainty Range
The uncertainty was evaluated by fitting the model parameters using i) a H 2 /D 2 dose increased or decreased by 30% and ii) the in-plane partition function and ZPE from either the present PBE PES or the RPBE PES from Ref. 17 for Pd(111). From Table S5, we estimate eV and ads = 1.02 ± 0.03 abs eV. = -0.093 ± 0.005

Temperature Dependent Rate Constants as Extended Arrhenius
To provide easy access to the isotope-specific rate constants, we parametrized the present model for recombination (Eq. S8 of the main paper) by fitting an extended Arrhenius equation (Eq. S28) to it in a temperature range from 300 to 1100 K. The extended Arrhenius fit shows typical deviations <1% with respect to the full model. The maximum deviation is 3%. The results are provided in Table S6. S17 Eq. S28 The equilibrium constant was parametrized with the same formula. Deviations do not A = A 01 / A bulk exceed 3%.

Reanalysis of Conrad et al. Isotherms
In order to quantify the recombination rate constants which underlie the work function isotherm data of Conrad et al. 19 an Arrhenius expression was used to simulate their data. In these experiments all coverages between 1 ML and 0 will be present depending on temperature and H 2

Eq. S30
Inclusion of bulk population does not affect the result and is not considered in this simulation. The resulting data of steady state hydrogen coverage at given temperature and pressure are compared to the experimental isotherms in Figure S7. A global amplitude parameter is used to transfer work function change to coverage. This assumes a linear relationship between work function and coverage. Figure S7: Fitting of adsorption Isotherms of Conrad et al. 19 The best fit parameters are , and . The small coverage a = 0.89 eV = 1.3 × 10 11 (MLs) -1 = 0.02 eV ML dependence is consistent with findings of Conrad et al. 19 12.

Reanalysis of Gdowski et al. TPD Data
In order to quantify the recombination rate constants which underlie the TPD data of Gdowski et al. 21 an Arrhenius expression was used to simulate their data. In TPD experiments all coverages between ~1 ML and 0 occur during the temperature ramp. Therefore, Eq. S29 was used here as well. The TPD data were simulated by propagating the following differential equation in time Eq. S31 The simulation was started at an initial coverage , which was obtained by the integral of the [H] 0 experimental traces, whereby the coverage of the 3.45 L dosage (see Figure 1 of Ref. 21 ) was set to 0.76 ML as the area of the TPD traces was not converged yet. The temperature at the beginning was set to 80 K and it increased with time according to the experimental heating rate The = 5.8 K s . S19 simulation leads to the time evolution of the hydrogen coverage and the TPD trace is obtained [H] as the H 2 production rate and is plotted as a function of temperature using .

= [H] 2 ( ) = 80K +
A global amplitude parameter is used to transfer the arbitrary experimental rate to the absolute simulated rate. The resulting fits are shown in Figure S8. Figure S8: Fit to TPD data of Gdowski et al. 21 The best fit parameters are , and . The resulting a = 0.91 eV = 1.

Tracer Kinetic Monte Carlo (TkMC) Method
To visualize H atom bulk penetration in the recombinative desorption experiment, we implemented a single particle kMC method. This allows to track the maximum penetration depth (MPD) of the propagated atom and correlate it to its time of recombination. This leads to the results presented in  The simulation proceeds as follows: For each trajectory, a is drawn from the experimental dosing 0 flux distribution. Then, the trajectory is propagated as follows: 1) Make a list of possible processes for the H atom in the current position. For example, at 0 these are recombination with the mean-field and hopping into the first subsurface.
2) Choose and execute a process according to principles of kMC 23 using random number .  Steps 1) to 4) are executed until the H atom undergoes recombination or the simulated time exceeds 2 ms. Figure S9 illustrates the propagation of a H atom. The simulation is carried out in a grid more coarsely grained than the one used for simulation of the mean-field kinetics to reduce simulation time.

Computational Details
The DFT calculations were performed with the Vienna ab-initio simulation package (VASP5.3.5). [24][25][26][27] For the description of electronic exchange and correlation effects, the PBE functional has been used for all calculations. 28,29 The optimum lattice constant of palladium was found to be 3.934 Å. The Pd(111) surface is modelled by a slab with 6 layers, whereas the two bottom layers were fixed.

× 2
Relaxation of all other layers was performed until forces and energies per atom were smaller than 1 meV/Å and 0.01 meV, respectively. To prevent interactions of the slab with its periodic images, a vacuum layer with a thickness of 14 Å was applied. The Brillouin zone was sampled by a 10 × 10 × 1 k-point grid making use of the sampling scheme proposed by Monkhorst and Pack. 30 For all calculations, the energy cut-off was set to 400 eV. The electronic energies of the H/Pd(111) PES were acquired by calculating 10 high symmetry points on the irreducible part of the surface. We permitted relaxation only of the H-atom's degree of freedom perpendicular to the surface plane (along theaxis) and the structure optimization was stopped when forces < 0.01 eV/Å were reached. The optimized geometry served as a starting point for the calculation of the harmonic H(D)-Pd stretch frequency, which was determined from the dynamical matrix using a finite difference method with four displacements along the -axis. The H(D)-Pd harmonic stretch frequencies were used for theaxis ZPE correction, see SI sec. 7.
The Pd(332) surface was modelled as a slab with four layers. The same convergence criteria (4 × 1) as for Pd(111) were applied for the relaxation of the slab. However, a k-point grid was used 4 × 4 × 1 to represent the Brillouin zone. The Pd(332) PES was sampled by calculating interaction energies, arranged in a 5×28 equidistant grid along the irreducible part of the Pd(332) surface. The corresponding harmonic frequencies were calculated in the same manner as for Pd(111).
In addition, the absorption energy was calculated using the PBE and the RPBE 31 functional. We represented the Pd bulk as super cell with being the optimal conventional lattice 2 0 × 2 0 × 2 0 0 constant of the respective functional (PBE: 3.934 Å, RPBE: 3.975 Å). In order to acquire the absorption energy, we placed an H atom at the interstitials of interest and subsequently let the super structure relax, where the cell volume was subject to optimization too. For integration over the Brillouin-zone, we utilized a Monkhorst-Pack k-point mesh. The energy cutoff of the 10 × 10 × 10 plane waves was set to 400 eV.